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Abstract 

The aim of this paper is to present a streamlined and fully three-dimensional version of 
the quasicontinuum (QC) theory of Tadmor et al. lITSlfTQl and to analyze its accuracy and 
convergence characteristics. Specifically, we assess the effect of the summation rules on 
accuracy; we determine the rate of convergence of the method in the presence of strong 
singularities, such as point loads; and we assess the effect of the refinement tolerance, 
which controls the rate at which new nodes are inserted in the model, on the development 
of dislocation microstructures. 

Keywords: A. Quasicontinuum; B. Error analysis; C. Nanoindentation; D. Atomistic mod- 
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1 Introduction 

The aim of this paper is to present a streamlined and fully three-dimensional version of 
the quasicontinuum (QC) theory of Tadmor et al. [>18iil9il and to analyze its accuracy and 
convergence characteristics. The theory of the quasicontinuum furnishes a computational 
scheme for seamlessly bridging the atomistic and continuum realms. The chief objective 
of the theory is to systematically coai^sen an atomistic description by-and only by-the 
judicious introduction of kinematic constraints. These kinematic constraints are selected 
and designed so as to preserve full atomistic resolution where required, e. g., in the vicinity 
of lattice defects, and to treat collectively large numbers of atoms in regions where the 
deformation field varies slowly on the scale of the lattice. Thus, in its purest form all input 
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into the theory concerning the behavior of the material is atomistic, and all approximations 
are strictly kinematic in nature. 

Different variants of the theory have been developed and documented over a series of 
publications U^ill[ll[l3|9l[ni[T3l[l7l[lll[l6l, where numerous examples of appli- 
cation have also been presented. Here we endeavor to emphasize the essential building 
blocks of the static theory, to wit: i) the constrained minimization of the atomistic energy 
of the solid; ii) the use of summation rules in order to compute the effective equilibrium 
equations; iii) and the use of adaption criteria in order to tailor the computational mesh to 
the structure of the deformation field. In particular, we develop a new class of summation 
rules based on sampling lattice functions over clusters of atoms. In Section |4j we present 
a detailed numerical analysis that probes various aspects of the accuracy and performance 
of the method. Specifically, we assess the effect of the cluster size and lumping proce- 
dure on accuracy; we determine the rate of convergence of the method in the presence 
of strong singularities, such as point loads; and we assess the effect of the refinement 
tolerance, which controls the rate at which new nodes are inserted in the model, on the 
development of dislocation microstructures. 

2 The quasicontinuum method 

We consider a reference configuration of a crystal such that its atoms occupy a subset 
of a simple d-dimensional Bravais lattice spanned by basis vectors {ai,i = 1, . . . , d}. 
The coordinates of the atoms in this reference configuration are: 

d 

X(Z) = ^rai, leCcZ'^. (1) 

i=l 

Here, I are the lattice coordinates which designate individual atoms, Z is the set of integer 
numbers and d is the dimension of space. From the definition, it follows that the set C 
is the collection of lattice sites occupied by atoms. The coordinates of the atoms in a 
deformed configuration of the crystal are denoted {q(/), I S £}. For convenience, we 
shall collect all atomic coordinates in an aiTay q and regard such an^ay as an element of 
the linear space X = R^'^, the 'configuration' space of the crystal. 

The energy of the crystal is assumed to be expressible as a function E{q), e. g., 
through the use of empirical interatomic potentials. In addition, the crystals may be sub- 
jected to applied loads. We assume these loads to be conservative and to derive from an 
external potential ^^^^{q). The total potential energy is, therefore: 

^{q)=E{q) + ^'^\q) (2) 

The crystal may also be subject to displacement boundary conditions over part of its 
boundary, e. g., as a result of the application of a rigid indentor. We wish to determine the 
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stable equilibrium configurations of the crystal. The problem is, therefore, to determine 
all the local minima of ^*(q) consistent with the displacement, or essential, boundary 
conditions. We shall, somewhat equivocally, enunciate this problem as: 

min$(g) (3) 

We emphasize, however, that in general the aim is not simply to determine the abso- 
lute minimizer of <&(g), but rather the much richer, and physically more relevant, set of 
metastable configurations of the crystal. 

The distinction arises from the fact that the energy function E{q) must be invariant 
under the class of affine deformations which map the crystal lattice onto itself. This im- 
plies, for instance, that E{q) must be periodic under crystallographic slip with period 
equal to the Burgers vector. Hence, the energy function E{q) is strongly nonconvex and 
the potential energy <&(q) in general has vast numbers of local minimizers, or metastable 
configurations. Techniques for systematically exploring the full set of metastable con- 
figurations, and indeed the entire configuration space of the crystal, include statistical 
mechanics and path-integral methods. An alternative approach is to load the crystal in- 
crementally and proceed by continuation. Thus, at every step in the solution procedure 
the loads, or the prescribed displacements, are incremented and the system is allowed to 
relax to a nearby stable configuration. The relaxation process may be governed by inertia, 
viscosity, kinetics, or some other rate-limiting mechanism; or may be numerical in nature 
and owe to the use of iterative solvers such as conjugate gradients. This latter approach is 
adopted in all the numerical tests presented in this paper. 

2.1 Interpolation 

The essence of the theory of Tadmor et al. lITSl [T9l is to replace Q by a constrained 
minimization of ^{q) over a suitably chosen subspace Xh of X. In order to define Xh, 
we begin by selecting a reduced set £/i C £ of A'^^j < 'representative atoms' (Fig.[l)- 
The selection of the representative atoms is based on the local variation of the fields and 
will be addressed subsequently. In addition, we introduce a triangulation of Ch. It 
bears emphasis that the triangulation Th may be unstructured. In particular, Ch need not 
define a Bravais lattice. The positions of the remaining atoms are determined by piecewise 
linear interpolation of the representative atom coordinates. We shall regard the resulting 
coordinates {qr^(Z), Z G £} as belonging to a linear space Xh of dimension Nhd. 

Let Lph{X\lh), Ih £ -C/i, be a collection of shape functions for Th- Thus, iph{X\lh) is 
continuous and piecewise linear, its domain is restricted to the simplices K ^Th incident 
to X{lh), and it vanishes at all nodes of the triangulation except at X{lh), where it takes 
the value 1, i. e., 

^h{X{l'h)\lh) = 5{l'h\lh). (4) 
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Figure 1: Example of triangulation Th of the crystal. 



By construction, 

Qhil) = Yl 'Ph{l\lh)qh{h), (5) 

where we write 

fh{l\h) = Vh{X{l)\h). (6) 

Evidently, {iph{l\lh)^ 1^ G Ch] constitutes a basis for and the fields q^il) are entirely 
determined by their values qh{lh) at the representative atoms. In addition, the basis lattice 
functions are required to satisfy the identity: 

i. e., the basis lattice functions must define a partition of unity over C. This requirement 
ensures that constant fields aie inteipolated exactly by the basis lattice functions. 

2.2 Reduced equations 

The reduced counterpart of problem ^ is now 



min ^{q^). (8) 
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The minimizers of the reduced problem satisfy the reduced equations of equilibrium 

fh{lh) = Y.f(^\lh)Ml\lh) = 0. (9) 

Here, 

/(g) = $,,(q) (10) 

are the forces corresponding to q and f{l\q) is the value of f{q) at site I. Thus, the 
reduced problem entails the solution of the N^d equations ^ in the N^d unknowns q{lh), 

2.3 Node-based summation rules 

The practicality of the method further hinges on the possibility of avoiding the calculation 
of the full atomistic force aiTay / and carrying out full lattice sums, as seemingly required 
in As noted by Tadmor et al. tl8 .191 . this may be accomplished by the introduction 
of summation rules similar to the conventional quadrature rules of numerical integration. 
The problem is thus to approximate sums of the general form 

S = Y,9{1), (11) 

where g{l) is a lattice function. By analogy to numerical quadrature rules, we begin by 
exploring summation rules of the form 

S^Y. ^h{l)9{l) ^ Sh (12) 

for some suitably chosen collection of summation points Sh and weights nh{l), not nec- 
essarily integer. Loosely speaking, the weights n/i(Z) may be regarded as the number of 
atoms represented by the site Z. Proceeding as in the development of numerical quadrature 
rules, the weights n/j(Z) may be determined by the requirement that the summation rule 
(fT2l be exact for a restricted class of lattice functions |13|. 

The lowest-order summation rule is obtained by requiring that all lattice functions in 
Xh, i. e., all piecewise linear- functions supported by the triangulation T^, be summed 
exactly. This is tantamount to setting Sh = Ch and requiring that the summation rule (IT2l 
be exact for all shape functions Lph{l\lh), Ih G ^h- This requirement gives, explicitly, 

nh{lh) = 'Yiph{l\lh), h G Ch- (13) 
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In fine regions of tlie triangulation, tiie sums on tlie rigiit-iiand side of (I13t may be com- 
puted explicitly. By contrast, in coarse regions of the mesh approaching the continuum 
limit this explicit calculation becomes impractical. However, in such regions the lattice 
sum (fT2l ostensibly reduces to an integral and the coiTcsponding weights are those of con- 
ventional Lobatto quadrature 13. In this limit, each simplex K ^ simply contributes 
{N/V)\K\/{d + 1) atoms to each of its d + 1 nodes, where N/V is the atom density in 
the undeformed configuration of the crystal. 

As a simple illustrative example, consider a monoatomic chain discretized into 'el- 
ements' each containing L bonds. For an internal node, the summation weight follows 
as 

n;, = 1 + 2 ( 1 - - J = L, (14) 

1=1 ^ ^ 

which coincides with the continuum limit. By contrast, for an end node the summation 
weight is 

- = '+|"(i-i) = ^- 

We note that Uh tends to the continuum Lobatto limit of L/2 as L oo and reduces to 
n/j = 1 in the full atomistic limit, as required. 

2.4 Reduced equations based on summation rules 

The application of a summation rule to the evaluation of the equilibrium equations (|9j 
leads to the system of equations 

fhih) « Yl Mi)fm)Mi\ih) = 0. (16) 

As desired, the calculation of the effective forces in (fT6t is of complexity 0{Nh) provided 
that the number of sampling sites in Sh is of order and the atomic interactions are 
short-ranged. 

The formulation of the method takes a somewhat unexpected turn at this point in that 
the node-based summation rule discussed in the foregoing suffers from a rank-deficiency 
problem. Thus, in the linear case the node-based summation rule leads to rank-deficient 
systems of equations. This is in analogy to the rank-deficiency of finite-element stiffness 
matrices computed by Lobatto quadrature, and manifests itself as the existence of so- 
called zero-energy modes which pollute the displacement field or render the system of 
equations singular altogether. 
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Figure 2: a) Cubic Lennard- Jones crystal used in numerical stability tests, b) 
Zero-energy deformation mode 



By way of illustration, we consider a Lennard-Jones crystal in the form of a cube 
containing 32 x 32 x 32 fee unit cells. In this example, eight representative atoms are 
introduced at the corners of the cube, and the cube is triangulated as shown in Fig.|2^. The 
stability of the crystal in its reference configuration may be assessed by examining the 
spectral properties of the reduced stiffness matrix. A direct calculation of its eigenvalues 
reveals that seven of them are identically equal to zero, which in turn denotes the presence 
of one zero-energy deformation mode. Fig. |2j). However, the zero-energy deformation 
mode is absent when the lattice sums ai^e canied out in full, i. e., when the test is repeated 
with the equilibrium equations (IT6t replaced by 

The numerical test just described shows that the node-based summation rules are in- 
deed rank-deficient, but it also suggests that the reduced system is stable provided that a 
sufficient number of sampling points is included in the summation rules. A new set of 
summation rules which offers this added flexibility is presented next. 

2.5 Cluster summation rules 

A class of summation rules which generalizes the node-based rule may be obtained by 
sampling the lattice function over neighborhoods of the representative atoms. We shall 
refer to these neighborhoods as clusters, and the resulting summation rules as cluster 
summation rules. Each cluster may be regarded as a representative crystallite where the 
state and the behavior of the crystal are sampled. 

More specifically, let C{lh) = {I : - X{lh)\ < r{lh)} be the cluster of lattice 

sites located within a sphere of radius r(Z/j) centered on the representative atom (Fig.|3ll- 
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Figure 3: Clusters of atoms in triangulation Th of the crystal. 



Assume for now that the clusters thus defined do not overlap. The cluster summation rule 
induced by the set of clusters C{lh),lh G is 

Sh= Mh)S{lh), (17) 

where S{lh) denotes the sum over all atoms in the cluster C{lh), i- e., 

S{h)= 9(1)- (18) 

iec{ih) 

As before, the cluster weights 7ih{lh),lh G ai'e computed by requiring that the cluster 
summation rule (fTTl be exact for all basis functions. 

It should be carefully noted that, as the triangulation size approaches the atomic length 
scale, the spherical clusters defined previously in general overlap, Fig.|4] In this case, we 
truncate the clusters as follows. A site I belongs to cluster C{lh) if its distance to X{lh) 
is less than a prescribed value r{lh) and less than its distance to any other representative 
atom. Thus, the cluster C{lh) is the intersection of the Voronoi cell containing X{lh) and 
sphere of radius r^l^) centered at X{lh). Ambiguous cases corresponding to sites which 
ai^e equidistant from two representative atoms are resolved randomly. In the limit of full 
atomistic resolution, each cluster contains exactly one atom and the coiTcsponding weight 
is 1. Thus, the cluster summation rules reduce seamlessly to the exact lattice sum when 
the fully atomistic limit is attained. 
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Figure 4: Triangulation % with two overlapping clusters. The dashed line demar- 
cates the boundary between the clusters. 
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Figure 5: A part of the ID monoatomic chain (h = 14 and r = 4). 



2.5.1 Summation-rule error analysis for monatomic chains 

The accuracy of cluster summation rules can be analyzed simply in the case of a one- 
dimensional monatomic chain. In this case, £ = Z is the set of all sites in the chain. 
Without loss of generality, we take the coordinate of a site / G £ in the undeformed con- 
figuration of the chain to be X(Z) = /. For simplicity, we consider a uniform triangulation 
of £ of size h. Thus, the coordinates of the representative atoms are X{lh) = Ihh, G Z, 
where h is the simplex size. We consider clusters of uniform radius r < h/2, Fig.|5l 

The summation weights follow from the requirement that the basis functions tph{l\lh) 
be summed exactly, with the result, 

nhih) = Ih G Z. (19) 

1 + zr 

Thus, for the infinite monatomic chain, nh{lh) is simply the ratio of the number of atoms 
in a simplex to the number of atoms in a cluster. 
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In order to obtain the summation error to leading order, we consider the sum of a 
quadratic polynomial supported on two adjacent simplices in the reduced chain, namely. 



^ \ otherwise 
The exact sum of g{l) over all sites of the chain is 

lec i=-h 

The calculation of the approximate sum Sh requires three local cluster sums Sh{—h), 
SiM and S,,(h) 

6 

l=-h 

5,(0, = g,. = r(i±lKl±^, 

S.(ft) = E ? = " + + (24) 

The use of the cluster summation rules leads to the approximate lattice sum 

Sh = rihiShi-h) + Sh{0) + Sh{h)] = ^Mi±^[2r2 + r(l - 3h) + 3h% (25) 

The summation error is | — 5| . An upper bound on the eiTor is obtained by setting r = 0. 
Then, for /i » 1 we have \Sh — S\r=o ~ (4/3)/i'^, which shows that the summation rule 
is of third order. 

The effect of the cluster size on the accuracy of the summation rule may be understood 
as follows. For simplicity, we consider the case of h ^ 1, whereupon SITli simplifies to 

S ~ ^h^ (26) 



and (|25ll to 
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where we write ^ = r/h. In the regime of ^ ^ 1, dZTt simpUfies further and the corre- 
sponding relative error behaves as 



\Sh - S\ 2 + r 

^ TT^ ^^^^ 

which is a decreasing function of r. Thus, as expected, increasing r from increases the 
accuracy of the summation rule. In the regime of r ^ 1 and r ~ /i, the relative error 
behaves as 

^^^~i(2e-i)(e-i) (29) 

which vanishes for = 1/2, or r = h/2. Thus, the exact sum is obtained when the 
clusters encompass the entire lattice, as required. 



2.6 Reduced equations based on cluster summation rules 



The application of the cluster based summation rule to the reduced equilibrium equations 
^ yields 



fhUh) 



rihil'h 



lecii'^ 



fil)Ml\h) 



(30) 



The computational complexity of (l30l is 0{NhNr), where A^^ is the number of lattice 
sites in a cluster of radius r. For A^^^. > 1, the computational effort is in excess of that 
corresponding to the node-based summation rules. Thus, the optimal value of the cluster 
size r{lh) is subject two opposing requirements. On one hand, the cluster size has to 
contain a sufficient number of atoms to ensure the appropriate accuracy and stability of 
the summation rule. On the other, however, such number should remain small for the 
method to retain its computational efficiency. The effect of the cluster size on the accuracy 
is investigated subsequently by way of numerical testing. 



2.7 Adaptive selection of representative atoms 

The third key component of the quasicontinuum method is the use of mesh adaption in 
order to tailor the computational mesh to the structure of the deformation field. Ideally, 
the adaptivity itself should be driven by the energetics of the system, i. e., the mesh should 
return the least possible potential energy for a fixed number of representative atoms. How- 
ever, in order to optimize the mesh in this manner it is necessary to know the relation 
between the energy, or suitable energy bounds thereof, and the mesh size. Unfortunately, 
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unlike the continuum case, an approximation theory for discrete systems such as consid- 
ered here appears to be entirely lacking at present, and such energy bounds are presently 
unavailable. 

Under these conditions, it becomes necessary to resort to empirical adaption indica- 
tors. In calculations we adopt as adaption indicator e{K) for simplex K the variation of 
the displacement field over K modulo rotations. This variation should be well-defined 
since the displacement field of the crystal may reasonably be expected to be of bounded 
variation. Specifically, we measure the variation of the displacement field as 



e{K) = y/\IlE{K)\h{K) (31) 

where IIe{K) denotes the second invariant of the Lagrangian strain tensor in simplex K, 
and h{K) is the size of K (cf IITsIITqI '). It follows from its definition that e{K) is invariant 
under rotations. The element K is deemed acceptable if 

< TOL (32) 



for some prescribed tolerance TOL < 1, and is targeted for refinement otherwise. In (I32t 
b denotes the magnitude of the smallest Burgers vector of the crystal. 

In our implementation of the method, the refinement of the simplices which violate 
d32t is accomplished by the application of subdivision rules. In particular, we choose 
to bisect the longest edge of the simplices which are targeted for refinement. The new 
representative atoms are inserted in the lattice site nearest to the geometrical midpoint of 
the longest edge. In order to preserve the quality of the computational mesh, we apply 
standard mesh-improvement operations such as edge-face or octahedral swapping |3JSQ 
and smoothing fTSll . Each new mesh is equilibrated and the remeshing criterion (l32t is re- 
evaluated for the new solution. The process is repeated until the mesh remains unchanged. 

The significance of d32b becomes apparent by considering a processes of crystallo- 
graphic slip across K. To this end, imagine that a pair of representative atoms in K 
undergo a relative sliding displacement of magnitude b across a slip plane of the crystal. 
For this deformation one has 

<K) = ^ (33) 

It follows, therefore, that the adaption criterion d32t ensures that the mesh size h{K) < b 
under the conditions just described. Thus, the adaption criterion is designed so that full 
atomistic resolution is attained when a simplex slips by a full Burgers vector. Interestingly, 
we have found that the dislocation patterns predicted in calculations are sensitive to the 
choice of tolerance, and generally TOL must be chosen much smaller than 1 in order not 
to inhibit dislocation nucleation. 
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3 Details of the computer implementation 



In order for the quasicontinuum method to perform optimally, some implementational 
issues require careful attention. Several potential performance bottlenecks are discussed 
in this section. 

3.1 Site-element mapping 

The numerical solution of the quasicontinuum equilibrium equations il6l . requires re- 
peated computation of atomic coordinates Qhil), I G J~- which depend, through the inter- 
polation constraints on the coordinates of representative atoms q^(Z/j), Ifi £ Ch. It 
is therefore important to devise an efficient algorithm for locating the simplex K £ Th 
which contains a given site I of the crystal. 

The simplex K £ containing site I £ C may be found from elementary geometry 
and an exhaustive search over 7/j. Such method carries, however, substantial computa- 
tional cost, due to a large number of required floating-point operations. Our approach is 
based on the fact that once the pair {I, K} has been established, it remains valid for the 
lifespan of 7/j. Thus, the following two stage strategy may be adopted: 

1. A general search routine, which for a given lattice site I £ C and triangulation 7/j 
returns K £ Tfi containing I. 

2. A look-up table (cache) that stores already associated pairs {/, K}. 

At first, the look-up table contains no {I, K} pairs and most of the inquiries results in 
the execution of the general search routine. As the force calculation proceeds, pairs are 
continuously inserted into the table, and fewer calls to the general routine are needed. At 
the end of the first force calculation, most of the pairs can be located in the cache, where 
they remain as long as Th stays unchanged. 

Our implementation of the look-up table is based on hashing, which is the method 
of referencing records in a table by doing arithmetic transformations of keys into table 
addresses. Any hashing algorithm requires two design decisions to be made: 

1. A hash function h{k) taking a key k as its argument and returning an index into a 
table must be chosen. 

2. One must establish a strategy for dealing with cases of two distinct keys ki and /c2 
for which the resulting values of the hash function are the same, i. e. h{ki) = h{k2) 
(collision resolution). 

In his monograph, Knuth provides an excellent description of hashing, including some 
possible choices of hash functions. For a simple three-dimensional Bravais lattice Q, the 
lattice coordinates of a site I = {/^, Z'^} £ may be considered as a key. The hash 
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function acting on this key produces an index into a table of simplices where the K, 
containing I, can be found. In cases in which a key contains more than one word (each 
of can be considered as a separate word) Knuth |6l suggests using the following hash 
function 

h{l) = [hi{l^) + h2{f) + h^{l^)] mod M. (34) 

where hi{k) denotes a hash function for and M is an integer parameter. The function 
a mod p computes the remainder of a /p. In general, the optimal choice of hi depends 
on Th and cannot easily be established. We have found, however, that the particular hash 
function 

/ii(r) = lshift(r, 2*-^), i = 1, 2, 3. (35) 

works well for a wide range of triangulations. The meaning of lshift(/,77T,) is simply 
"shift all bits in I left m positions". The operations involved in the calculation of hi are 
elementary, and the computational cost of each such calculation is small. The integer 
parameter M in (l34l i may be chosen arbitrarily, but Knuth ||6l points out that, when it 
is taken to be some power of 2, i. e., M = 7f\ a mod M is equivalent to masking the 
low p-bits from a (on most of cuiTcntly available computer hardware). The application 
of (l34t to a lattice site I results in an index into the hash table from the interval [0, 2^). 
Accordingly, the size of the table is limited to 2^, which renders collisions likely. The 
collisions are handled with the help of additional short tables which store {l,K} pairs 
for all I such that the h{l) produces the same value. Therefore, the search through the 
cache for a simplex K containing site I becomes essentially a two stage process, in which 
the application of the hash function is followed by a linear search. Since, the number of 
lattice sites having the same index is small, and the additional cost due to the linear search 
phase is not significant. 

3.2 Computation of cluster weights 

As explained earlier, the cluster weights nh{lh), I ^ Ch obtained requiring that the 
shape functions are summed exactly by the cluster summation rule dTTl . The applica- 
tion of the summation rule to shape functions (ph{l\lh) leads to the system of Nh linear 
algebraic equations 

^ AihlQ n(Z'J = b{h), h G Ch, (36) 

to be solved for nh{lh), I G C-h- Here the matrix A{lh\l'i^), I'h^lh G is 

Ailh\l'h)= E Ml\h)- (37) 

Z'GC(i;) 
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and C(Zjj) denotes the cluster of lattice sites centered at l'^. The array b{lh), Ih is simply 

b{lh) =Y,Ml\h)- (38) 

The calculation of the aiTay b requires the evaluation of the shape functions at all lattice 
sites within a crystal. While this may be regarded as acceptable for small crystals, it 
becomes prohibitively expensive when appUed to large samples. However, as we have al- 
ready shown, in the continuum limit summation may be replaced by integration, in which 
case simplex contributes N{K)/{d+l) to each of its vertices. N{K) = {N/V)\K\ is 
the approximate number of sites within K and N/V is the atom density in the undeformed 
configuration of the crystal. Specifically, one may introduce a cutoff Nc and restrict the 
direct calculation of the sum in d38b to simplices K for which N{K) < Nc- 

It is apparent from (l37l that A has a sparse structure, which suggests the use of spe- 
cialized solvers for sparse linear systems. However, a simpler alternative route is to resort 
to lumping in order to replace A by a diagonal matrix (see, e. g., Q). A lumping tech- 
nique which is widely used in finite elements is the row-sum technique, which gives the 
diagonal entries of the lumped matrix as 

A{h\h)= E '^(^l^'^)' heCh, (39) 

with all other entries in A set to zero. Once the matrix A is lumped, the solution of (l37l 
is trivial and gives 

"(^'^) = TTTTTT' ^ ^f^- (40) 

A{lh\lh) 

In the atomistic limit, the values of cluster weights computed from (l40l become identical 
with obtained exactly. 



4 Numerical tests 

The accuracy of the quasicontinuum method is largely determined by three factors: 

• The value of the cluster-size parameter r, which controls the accuracy of the cluster 
summation rules. 

• The approximations introduced in the computation of the cluster weights, namely: 
the approximate calculation of vector b, controlled by the cuttoff parameter Nc, 
and the lumping procedure for constructing the diagonal matrix A. 
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Figure 6: Crystallographic orientation of test sample used in the simulations of 
nanoindentation. 



• The value of parameter TOL in (l32l . which controls the process of representative 
atom insertion. 

In this section we present an analysis of the influence of each of these factors on the 
accuracy of the quasicontinuum method. We measure the quality of the solution in energy 
terms. Specifically, we identify the eiTor in the approximate solution with 



where q is the solution of the full system obtained, e. g., by relaxing q^. Since the 
potential energy decreases during this relaxation, the energy error is greater or equal to 
zero. For a well-supported harmonic crystal the energy en^or defined in dTO defines a 
proper norm of the eiTor lattice function q^ — Q- 

4.1 Test problem definition 

We take nanoindentation as a convenient test problem for assessing the performance of 
the method. A salient feature of nanoindentation is the presence of a highly nonuniform 
state of deformation resulting in a sharp mesh gradation away from the indentor This 
feature effectively tests the adaptivity of the method. 



(41) 
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Figure 7: Cross section through the quasicontinuum test sample used in the 
nanoindentation simulations. A part of the indentor is also shown. The white 
circle represents a cluster of radius r = -\/69[o"] (corresponding to sixty-four shell 
of neighbors). 



17 



The test sample used in simulations is a cube of an fee nearest-neighbor Lennard- 
Jones crystal containing 64 x 64 x 64 fee unit cells, or a total of 1 073 345 atoms. The size 
of the cube is limited by the need to compute the solution of the full problem by direct 
atomistic simulation in order to evaluate the eiTor. The Lennard- Jones potential is Cl|8l| 



where a and e are parameters which set the length and energy scales, respectively. The 
choice of the Lennard- Jones potential is motivated by the desire to eliminate surface ef- 
fects which would otherwise compound the interpretation of results. 

The surfaces of the sample are aligned with the cube directions. Fig. |6l Atoms in 
the bottom surface are constrained to remain at their initial positions throughout the test, 
whereas in all side surfaces the atoms are allowed to move in the z-direction only. No 
displacement constraints are introduced on the top surface of the cube. 

In calculations we use a model of a spherical indentor proposed by Kelchner et al. 13. 
In this model, the indentor is regarded as an additional external potential interacting 
with atoms in the substrate. The potential has the particular functional form 



where R is the radius of the indentor, r denotes distance between a site and the center of 
the indentor, ^4 is a force constant and H{r) is the step function. In calculations we adopt 
the following values of the parameters: R = 100 [a] and A = 2000 [e/cr^]. 

The initial triangulation of the cube is specifically tailored to the nanoindentation ge- 
ometry, Fig.0 Thus, in the region of the crystal located directly underneath the indentor, 
full atomistic resolution is introduced from the outset. Away from this region, the trian- 
gulation becomes gradually coarser. The resulting number of representative atoms in the 
initial mesh is 888, a significant reduction from the total number of atoms (1 073 345) in 
the sample. A cross section of the initial mesh through the center of the cube with a plane 
y = const is show in Fig.El AH solutions ai^e computed using a nonlinear- version of the 
Conjugate Gradient method. 

4.2 Effect of the cluster size 

The cluster radius r controls the accuracy of the cluster summation rules. In order to 
isolate the effect of the cluster size from other factors, in the present tests the summation 
weights are computed exactly and the mesh is kept fixed. The indentor is pushed through 
a distance 0.1 [a] into the crystal along the z-axis. This value ensures that the material 
immediately under the indentor becomes highly deformed and is pushed well into the 
nonlinear regime. However, the induced deformation is not sufficient for any defects to 




(42) 



$^^t(r) = AH{R-r){R-rf 



(43) 
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Figure 8: The error e as a function of the cluster size r. 



arise. This greatly simplifies the interpretation of results and the error analysis, since 
under the stated conditions the full atomistic solution q to compare with is well-defined. 

The computed dependence of the error e on the cluster size r is shown in Fig. [8l It 
is evident from this plot that the inclusion of a single shell of neighbors in the clusters 
suffices to eliminate the rank-deficiency of the node-based summation rule. As suggested 
by our previous analysis, the en^or is greatest for the smallest cluster size and decreases 
rapidly as the cluster size increases. For r above \/8 [a] , con^esponding to the radius of the 
eighth shell of neighbors, the solution becomes insensitive to the cluster size. Thus, the 
use of relative small clusters results in high accuracy while preserving the computational 
efficiency of the method. 



4.3 Effect of the lumping procedure 

The enor in the calculation of the cluster weights originates from two sources: the approx- 
imate computation of the array b in (l38l . controlled by the value of parameter N^, and the 
lumping of matrix A in (I37t . In order to appraise this error, we repeat the calculations 
described in the preceding section with the summation weights computed approximately. 
The cutoff Nc is set to 2000, i. e. the contributions to the aiTay b from simplices which 
contain fewer than 2000 atoms are computed explicitly. The dependence of the error 
e on the cluster size parameter r for both the approximate (squares) and exact (circles) 
weights is shown in Fig. |9l Small differences do arise for very small clusters, and virtu- 
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Figure 9: The error e as a function of the cluster size r. Squares and circles corre- 
spond to the approximate and exact calculation of cluster weights, respectively. 



ally disappear for large clusters. In view of these results we may conclude that the use of 
approximate summation weights does not result in a significant loss of accui^acy. 

4.4 Convergence properties of the quasicontinuum method 

We proceed to assess numerically the rate at which the quasicontinuum solution converges 
towards the fully atomistic solution. The test case is as in the preceding sections. In all 
cases, the cluster size r is set to 2^/2 [a] and the cluster summation weights are com- 
puted approximately with a cutoff parameter Nc = 2000. The coarsest mesh is shown 
in Fig.fTOh. Subsequent finer meshes aie constructed by regular refinement, Fig.fTOb. re- 
sulting in increasing numbers of representative atoms A'^^, with the exception of the fully 
resolved case, h = \f2 \a\ , in which the representative atoms are placed on all the crystal 
lattice sites. In these calculations, the crystal is loaded simply by imparting a downward 
displacement —0.1 [a] to the central atom on the surface labeled A in Fig.|6l 

Fig. ^2 shows the variation of the energy error with the number of representative 
atoms. Three distinct regimes are evident in this figure. For small values of Nfi, the 
eiTor behaves as e ~ 0{Nh~°'), with a convergence rate a « 0.39. At Nh = 4913, 
conxsponding to h = 4\/2 [a], the summation clusters begin to overlap and the rate of 
convergence increases to a « 0.55. The final drop of the error occurs when full atomistic 
resolution is attained. 
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Figure 10: Examples of triangulations used in the convergence study. Plots corre- 
spond toh = 64^2 [a] (left) and h = 16^2 [a] (right). 




Figure 1 1 : The energy error e as a function of number of representative atoms N^. 
in the sample. 
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These results demonstrate the convergence characteristics of the quasicontinuum method 
in the extreme case of an applied point load. When full lattice sums are performed the 
calculated convergence rate of 0.55 approaches the optimal rate of convergence of linear 
finite elements applied to linear elasticity, which is 2/3. The introduction of cluster sum- 
mation rules results in a degradation of the convergence rate. It should be carefully noted, 
however, that some of these conclusions may well depend on the loading geometry. As is 
well known, the linear elasticity solution corresponding to a point load diverges under the 
point of application of the load and does not possess finite energy. The finiteness of the 
energy of atomistic solution owes entirely to the discreteness of the atomic lattice. In this 
sense, the test case considered here is, therefore, a worse case, and it is possible that the 
converge of the quasicontinuum method is more robust for smooth loading. 

4.5 Effect of the remeshing-indicator tolerance 

The lack of convexity of the energy functional for crystalline materials allows for lattice 
defects to develop. Among these defects, dislocations play a crucial role, as the carriers 
of the plastic deformation. One of the appeahng features of the quasicontinuum method 
is its ability to follow the nucleation and motion of dislocations inside the crystal. This is 
achieved by adaptively providing full atomistic resolution in the highly deformed regions 
of the crystal. In the nanoindentation test, dislocations nucleate under the indentor upon 
the attainment of a critical load and subsequently propagate into the crystal. 

The insertion of new representative atoms, leading to mesh refinement, is controlled 
by criterion d32t and, specifically, by the tolerance TOL. A small value of TOL promotes 
refinement, whereas a large value of TOL inhibits refinement. Most importantly, an ex- 
cessively large value of TOL may have the undesirable effect of inhibiting the nucleation 
of dislocations under the indentor altogether. On the other hand, an unduly small value of 
TOL results in excessive refinement and a prohibitive computational burden. 

In this section we explore this trade-off by way of numerical testing, with a view to 
bracketing the optimal value of the tolerance parameter TOL. To this end, we fix the 
cluster size to r = \/2 [a] , and we use the approximate cluster summation weights with 
a cutoff Nc = 2000. The indentor is driven into the crystal at increments of 0.1 [a] 
up to a maximum indentation depth of 2.0 [a]. The first loading step does not result in 
the nucleation of dislocations, even in the fully atomistic model, and is identical to the 
calculations described in the preceding sections. By contrast, at the indentation depth of 
2.0 [a] the presence of dislocation structures under the indentor is clearly evident in the 
fully atomistic model. Fig. 

The energy enw e is plotted in Fig.^Jas a function of indentation depth 5 for different 
values of the tolerance parameter TOL. Initially, the en^or is insensitive to TOL, and all 
curves ostensibly follow the same linear growth pattern. However, the error curves diverge 
after a penetration depth 5 = 0.5 [a]. For large values of TOL, the error continues to 
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Figure 12: Dislocation structure at the indentor penetration S = 2.0 [a] predicted 
by the fully atomistic model. The figure displays the energetic atoms (red) under- 
neath the crystal surface (blue). 
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Figure 13: Error e vs indentation depth 6 for different values of the remeshing 
tolerance parameter TOL. 



grow owing to the inability of the model to resolve the emerging dislocation structures. 
Indeed, the nucleation of dislocations during the test is entirely inhibited by all but the 
smallest values of the tolerance, TOL = 0.004. For small values of the tolerance, TOL = 
0.005 and 0.004, the energy eiTor reaches a plateau at roughly 1% of the total energy, and 
remains constant or decreases thereafter. 

The tinal mesh at 5 = 2.0 [a] corresponding to a tolerance of TOL = 0.004 is 
shown in Fig. El The total number of representative atoms in this configuration is = 
141, 469, which con^esponds approximately to 1/8 of the total number of atoms in the 
sample. It should be stressed, however, that since all of new representative atoms are 
inserted in the vicinity of the indentor, the ratio Nh/N can be made arbitrarily small by 
considering increasingly large samples. 

The dislocation structures predicted by the quasicontinuum method for a tolerance 
TOL = 0.004, Fig. should be contrasted with those predicted by the fully atomistic 
model. Fig. El The figures display the energetic atoms (red) underneath the surface of 
the crystal (blue). Evidently, these structures differ in detail. This is expected, owing to 
the massive lack of uniqueness which characterizes the problem. Thus, as dislocation are 
nucleated and the lack of convexity of the energy comes into play, numerous deforma- 
tion paths become available to the crystal. Many of these deformation paths differ only 
slightly in their energy content and are, therefore, indistinguishable at the macroscale. 
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Figure 14: The cross section through the quasicontinuum sample at the indentor 
penetration 6 = 2.0 [6]. 
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Figure 15: Dislocation structure predicted by the quasicontinuum model at inden- 
tation depth 6 = 2.0 [a]. The figure displays the energetic atoms (red) underneath 
the crystal surface (blue). 
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Figure 16: Load displacement curve predicted by atomistic (LS) and quasicon- 
tinuum (QC) simulations. 



e. g., in terms of the con-esponding force-depth of penetration curves. This is specially 
so when the crystal is loaded along a direction of high symmetry, which results in the 
possible activation of a large number of competing slip systems. Under these conditions, 
a pointwise comparison of solutions becomes mathematically meaningless, and the sole 
meaningful criterion to measure the quality of a solution is its energy contents. By this 
criterion, the quasicontinuum and fully atomistic solutions are of compai'able quality, as 
their energies differ by less than 1%. 

This contention is confirmed by Fig. which compares the force on the indentor 
/ as a function of penetration depth 5 predicted by the fully atomistic and quasicontin- 
uum (TOL = 0.004) models. Indeed, the two curves are ostensibly indistinguishable 
even beyond the onset of dislocation nucleation. On this basis, in practice choices of the 
tolerance TOL in the range of 0.001-0.004 would appear to strike an adequate balance 
between accuracy and performance demands. 



5 Summary and Discussion 

We have developed a streamlined and fully three-dimensional version of the quasicontin- 
uum method of Tadmor et al. fTSl [T9ll and we have presented a numerical analysis of 
its accuracy and convergence characteristics. As a new addition to the theory, we have 
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formulated a new class of summation rules in which the lattice function being summed is 
sampled over clusters of atoms. The size of these clusters is an adjustable parameter and 
controls the accuracy of the summation rule. We have also presented an efficient method 
for computing the requisite summation weights in linear time. Beyond its usefulness as a 
numerical scheme for approximating lattice sums, the cluster approach provides a system- 
atic means of sampling the behavior of small representative crystallites, and thus opens 
a possible avenue for incorporating additional physics to quasicontinuum models such as 
diffusion of solute atoms and thermal lattice vibrations. 

We have presented a suite of numerical tests which demonstrate the accuracy and per- 
formance of the method. As expected, the accuracy of cluster summation rules increases 
with cluster size. Furthermore, our numerical tests suggest that the addition of a sin- 
gle shell of neighbors suffices to stabilize the rank-deficiency which afflicts node-based 
summation rules. The computed convergence rate of the method in problems involving 
the appUcation of point loads to crystals is close to linear when the lattice sums are per- 
formed exactly, and decreases somewhat when the sums are approximated using a cluster 
summation rule. It is worth noting that, contrary to continuum problems for which a 
well-developed approximation theory exists, a similar approximation theory for lattice 
problems appeai^s to be missing at present, even for linear problems. The development 
of rigorous error bounds for finite-element approximations to lattice problems is a clear 
worthwhile direction of future research. 

Perhaps the most interesting issue among those addressed in this paper concerns the 
ability of the quasicontinuum method to simulate microstructural evolution, which owes 
largely to mesh adaptivity. The same lack of convexity which allows for defects and mi- 
crostructures to arise in the first place renders solutions massively nonunique. For a given 
loading or prescribed deformation, it is in general possible to find multiple equilibrium 
states of a crystal with defects possessing equal or nearly equal energies. These states ai^e 
indistinguishable from a macroscopic point of view, e. g., in the sense of yielding iden- 
tical force-displacement curves. Whether the solution follows one deformation path or 
another depends sensitively on small perturbations of the system, including details of the 
mesh design. Under these conditions, the pointwise comparison of solutions is not partic- 
ularly meaningful. Our numerical tests suggest that, with sufficient mesh adaptivity, the 
quasicontinuum method is capable of simulating evolving microstructures comparable, in 
energetic terms, to those obtained from a full atomistic calculation. 
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